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Using classical molecular dynamics simulation, we have studied the effect of edge-passivation by hydrogen (El- 
passivation) and isotope mixture (with random or supperlattice distributions) on the thermal conductivity 
of rectangular graphene nanoribbons (GNRs) (of several nanometers in size). We find that the thermal 
conductivity is considerably reduced by the edge H-passivation. We also find that the isotope mixing can 
reduce the thermal conductivities, with the supperlattice distribution giving rise to more reduction than the 
random distribution. These results can be useful in nanoscale engineering of thermal transport and heat 
management using GNRs. 



Graphene^] is a monolayer of graphite with a honey- 
comb lattice structure. It exhibits many unique proper- 
ties and has drawn intense attention in the past few years. 
The unusual electronic properties of graphene are promis- 
ing in many fundamental studies and applications, e.g., 
the ultrahigh electron mobility^ and the tunable band 
gap and magnetic properties by the size and edge chiral- 
ity of GNRsP^l Graphene also has remarkable thermal 
properties. The measured value of thermal conductivity 
of graphene reaches as high as several thousand of W/m- 
KpESl amon g th e highest values of known materials. Pre- 
vious studies^"^ show that the thermal transport in 
GNRs depends on the edge chirality of GNRs. In realistic 
graphene samples, the edges are often passivated^^ and 
the isotope composition can be controlled.^ Motivated 
by these, we study the effect of the edge H-passivation 
and various isotope distributions on the thermal trans- 
port in GNRs. We find that the thermal conductivity can 
be reduced by the edge H-passivation and tuned by the 
isotope distributions. Our study is useful in nanoscale 
control and management of thermal transport by engi- 
neering the chemical composition of GNRs. 

In this work, we employ the classical molecular dy- 
namics (MD, similar to the method in Ref. ITTj) to cal- 
culate the thermal conductivities of GNRs. We use the 
Brenner potential,^ which incorporates the many-body 
carbon-carbon and carbon-hydrogen interactions by in- 
troducing a fractional number of covalent bonds. This 
method has b een successfully applied to m any car bon- 
based systems} 19 * 20 ! especially to graphene! 11 * 21 * 22 ! The 
structures of GNRs are shown in Fig. [I] (with edge H- 
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passivation) and the insets of Fig. [3] (with isotope mix- 
tures). We use fixed boundary conditions, i.e., the atoms 
denoted by squares are fixed at their equilibrium posi- 
tions. The atoms denoted by left- a nd rig ht-pointing tri- 
angles are placed in two Nose-Hoove j 23 * 24 ! thermostats at 
different temperatures. The thermal conductivity can be 
calculated from the Fourier law n = Jdj \/STwh), where 
AT is the temperature difference (chosen to be in the 
linear response regime^) between two thermostats, J is 
the resulting thermal current, d is the length, w is the 
width and /i(=0.335 nm) is the thickness of GNRs, re- 
spectively. Calculations presented below are performed 
for representative GNRs with d ~ 6nm and w ~ 1.5 
nm. All the calculated thermal conductivities are nor- 
malized by kq. kq is the thermal conductivity calculated 
at 100 K for the pure 12 C GNR with armchair edge and 
without H-passivation (shown in Figure (a) in Ref. 
Although the specific value of k depends on the GNR 
sizepSl the choice of the rmostat and boundary condition 
in MD simulation,* 12 * 27 ^ we have checked that our conclu- 
sions and the qualitative behavior of k discussed below 
do not. The equations of motion for atoms labeled by 
triangles in either left or right Nose-Hoover thermostat 
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where i runs over all the atoms in the thermostat, p, is 
the momentum of the z-th atom, F.; is the total force act- 
ing on the z-th atom, 7 is the dynamic parameter of the 
thermostat with initial value of zero, r is the relaxation 

time of the thermostat, T(t) = 3P ? kB Yli ^7 i s tne actual 
temperature of atoms in the thermostat at time instant t 
, T = Tl (or Tji ) is the desired temperature of the left 
(or right) thermostat, N is the number of atoms in the 
thermostat, kg is the Boltzmann constant and rrii is the 
mass of the i-th atom. The atoms labeled by circles obey 
the Newton's law of motion with 7 = 0. The carbon- 
carbon potential is the same for all carbon isotopes. 

First, we study the temperature dependent thermal 
conductivity of H-passivated GNRs. Fig. [l] shows the 
armchair and zigzag GNRs with top and bottom edges 
H-passivated. As shown in Fig. [2] the edge H-passivation 
significantly reduces the thermal conductivity, compared 



FIG. 1. Structure of hydrogen-passivated armchair (a) and 
zigzag (b) GNRs. The hydrogen atoms are denoted by smaller 
symbols while the 12 C atoms are denoted by larger ones. ■ 
denotes fixed boundary atoms. ► (■<) denotes atoms in the 
left (right) thermostat. • denotes the remain atoms in the 
bulk. 
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FIG. 2. Temperature dependent thermal conductivity of 
GNRs with and without edge H-passivation. 

to the non-passivated GNRs. A recent studj@ using 
equilibrium MD has obtained similar conclusions. We 
note that the error bars (related to molecular dynam- 
ics fluctuations) for H-passivated GNRs are considerably 
larger than that for the non-passivated GNRs, probably 
due to the much smaller mass of hydrogen atoms. 

We also study the effect of the mixture of carbon 
isotopes 12 C and 13 C on the thermal conductivity of 
GNRs. Here, we demonstrate the results in the case 
of armchair GNRs (qualitatively similar results are ob- 
tained for zigzag GNRs). The concentration of 13 C is 
Ni 3 /(N 12 + N 13 ), where N 12/13 is the number of 12 / 13 C 
atoms. The thermal conductivity is seen to be reduced 
by introducing 13 C, and the thermal conductivity of pure 
13 C GNRs is lower than that of pure 12 C GNRs because 
13 C atoms have larger mass and thus give lower phonon 
frequencyPS The inset of Fig. shows a typical ran- 
dom isotope distribution. Another isotope distribution 
pattern we study is the isotopic supperlattice structure 
shown in the inset of Fig. [3jb) . Here, the whole GNR 
is composed of four slices with equal length and the 
same isotope composition. Within each slice (such as 
the dashed box in the inset of Fig. (3^b)), the number 
(L) of vertical 13 C atomic chains with zigzag shape (e.g., 
L = 4 for the inset of Fig. [3^b) ) can vary from to 7 
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FIG. 3. Temperature dependent thermal conductivity of 
GNRs with 13 C isotopes distributed (a) randomly and (b) in 
a supperlattice structure. The insets show the corresponding 
typical structures of GNRs with the same meaning of sym- 
bols as that in Fig. [T] The larger (smaller) symbols denote 
13 C ( 12 C) atoms. 



(see details in Ref. 26). L = (L = 7) corresponds to 
the pure 12 C ( 13 C) GNR. The temperature dependent 
thermal conductivity in Fig. [3] shows that the isotope ef- 
fect becomes more evident at higher temperatures. We 
show the thermal conductivity as a function of the con- 
centration of 13 C calculated at the temperature of 500 
K in Fig. |4| In the case of the random distributions, 
the calculated thermal conductivity is an average of 10 
different distributions with the same isotope concentra- 
tion. For random isotope distributions, the isotope con- 
centration dependent thermal conductivity (dashed line 
in Fig. [4]) shows a pan shape and is relatively flat in the 
concentration range of ^20-90%. The thermal conduc- 
tivity is reduced by ~10% around the isotope concen- 
tration of ^50%. The error bars are determined from 
the deviations of the thermal conductivities for the 10 
different distributions from their average value. In con- 
trast, the conical shape of the solid line in Fig. [4] shows 
much stronger dependence of the thermal conductivity 
on the isotope concentration for the supperlattice struc- 
tures, with ^30% reduction of the thermal conductivity 
at ^50% of the isotope concentration. We have also ob- 
tained similar result^Sl using velocity exchange MESS in 
the LAMMPS package P 

It has been previously demonstrated that the thermal 
conductivity depends on the edge chiralityS^K^ in the ab- 
sence of H-passivation, i.e., the thermal conductivity of 
the zigzag GNR is larger than that of the armchair GNR. 
However, as shown in Fig.[2j their thermal conductivities 
become close to each other within the MD error bars af- 
ter the H-passivation, suggesting that phonon scattering 
with the hydrogenated edges dominates over the contri- 
bution from the chirality effect. We have suggested that 
the smaller thermal conductivity of armchair GNRs is 
due to the stronger phonon scattering at the armchair 
edgesP^ It is interesting to note that the H-passivated 
zigzag GNR in Fig. [ljb) resembles the armchair GNRs 
at the edges. We suggest that the thermal transport in 
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FIG. 4. Thermal conductivity as a function of the 13 C con- 
centration for supperlattice (solid line) and random (dashed 
line) isotope distributions. 



small GNRs (several nanometers in size in our study) is 
strongly affected by the edge configuration. 

Recently, it has been experimentally demonstrated 
that different carbon isotopes can be controllably intro- 
duced in graphene, such as 13 C, in the chemical vapor 
deposition growth of graphene on metals. Both random 
and segregation (by domains of different isotopes) distri- 
butions have been observed.^ This opens possibilities of 
engineering the thermal properties of graphene by isotope 
distributions. The isotope effect on the thermal trans- 
port has been studied in several nanomaterial s, suc h as 
carbon nanotubesf^ boron nitride nanotubesj 34 * 35 ! and 
silicon nanowiresP^ The pan shape of the dashed line in 
Fig. [4] is consistent with the the reduction of the thermal 
conductivity in the "alloy limit. A similar pan shape 
is found in GNRs^ and SiGe nanowires^ by tuning the 
composition. By keeping the isotope concentration a con- 
stant of ~50%, it has been showrF^that the thermal con- 
ductivity as a function of the slice length (which is kept 
constant in our simulations) gives similar conical shaped 
dependence as we see in Fig. [4] for the supperlattice struc- 
tures. 

In conclusion, the classical MD is applied to calculate 
the thermal conductivities of rectangular GNRs. We find 
that the edge H-passivation can reduce the thermal con- 
ductivity significantly. We also show that the thermal 
conductivity depends on the concentrations of isotopic 
atoms and their distribution patterns. The isotopic sup- 
perlattice distributions can reduce the thermal conduc- 
tivity much more than random distributions. These find- 
ings can be useful in controlling heat transfer in nanoscale 
using GNR-based thermal devices. 
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SUPPLEMENTARY MATERIAL 

Supperlattice structures of different 13 C concentrations 
are shown in the following plots. The larger (smaller) 
symbols denote the 13 C ( 12 C) atoms. The number L, 
denoting the number of columns of zigzag 13 C atomic 
chains in each subplot, varies from to 6 for these plots. 
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(b) 13.1% 13 C, L=1 



0.! 



-0.! 




1 2 3 4 5nm 



(c) 27.8% 13 C, L=2 
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(e) 57.4% 13 C, L=4 
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(f) 72.2% 13 C, L=5 
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(g) 86.9% 13 C, L=6 
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